num.stems <- c(6,8,9,6,6,2,5,3,1,4)
sum(num.stems)
# data frame of tabulated data
aphids <- data.frame(aphids=0:9, counts=num.stems)
aphids
names(num.stems) <- 0:9
barplot(num.stems)
aphid.data <- rep(0:9,num.stems)
aphid.data
poisson.LL <- function(lambda) sum(log(dpois(aphid.data, lambda)))
poisson.negloglik <- function(lambda) -poisson.LL(lambda)
out.pois <- nlm(poisson.negloglik,3)
out.pois
 
# Poisson regression with log link
out.p <- glm(aphid.data~1, family=poisson)
summary(out.p)
coef(out.p)
exp(coef(out.p))
 
# Poisson regression with identity link
out.p1 <- glm(aphid.data~1, family=poisson(link=identity))
summary(out.p1)
 
# Poisson probabilities
dpois(0:9, lambda=out.pois$estimate)
# expected counts
dpois(0:9, lambda=out.pois$estimate)*50
# missing tail probability
sum(dpois(0:9, lambda=out.pois$estimate))
 
# include tail probability
c(dpois(0:9, lambda=out.pois$estimate), 1-ppois(9,out.pois$estimate))
sum(c(dpois(0:9, lambda=out.pois$estimate), 1-ppois(9,out.pois$estimate)))
 
# expected counts
c(dpois(0:9, lambda=out.pois$estimate), 1-ppois(9,out.pois$estimate))*50
 
# counts with extra category for X > 9
rbind(c(num.stems,0), c(dpois(0:9, lambda=out.pois$estimate), 1-ppois(9,out.pois$estimate))*50)
 
# counts in which last category is 9 or more
rbind(num.stems, c(dpois(0:8, lambda=out.pois$estimate), 1-ppois(8,out.pois$estimate))*50)
 
# expected Poisson counts
exp.pois <- c(dpois(0:8, lambda=out.pois$estimate), 1-ppois(8,out.pois$estimate))*50
exp.pois
 
# two different ways to calculate the last category
1-ppois(8, out.pois$estimate)
ppois(8, out.pois$estimate, lower.tail=F)
 
# bar plot with pooled last category
out.bar <- barplot(num.stems, names.arg=c(0:8,'9+'), ylim=c(0,11))
# x-coordinates of midpoint of bar
out.bar
# add expected counts to bar plot
points(out.bar, exp.pois, pch=16, type='o', cex=.9)
 
# expected probabilities
expected.p <- c(dpois(0:8, lambda=out.pois$estimate), 1-ppois(8,out.pois$estimate))
expected.p
 
# pool categories so all expected counts exceed 5
expected2 <- c(sum(exp.pois[1:2]), exp.pois[3:6], sum(exp.pois[7:10]))
expected2
# repeat process with observed counts
observed2 <- c(sum(num.stems[1:2]), num.stems[3:6], sum(num.stems[7:10]))
observed2
 
# Pearson statistic
sum((observed2-expected2)^2/expected2)
length(observed2)
# critical value of test statistic
qchisq(.95,6-1-1)
# p-value of test statistic
1-pchisq(sum((observed2-expected2)^2/expected2), 6-1-1)
 
# repeat using chisq.test function
exp.p <- c(dpois(0:8, lambda=out.pois$estimate), 1-ppois(8,out.pois$estimate))
exp.p
expected2.p <- c(sum(exp.p[1:2]), exp.p[3:6], sum(exp.p[7:10]))
expected2.p
chisq.test(observed2, p=expected2.p)
 
# without pooling chisq.test issues a warning
actual.p <- c(dpois(0:8, lambda=out.pois$estimate), 1-ppois(8,out.pois$estimate))
actual.p
chisq.test(num.stems, p=actual.p)
 
# parametric bootstrap
out.sim <- chisq.test(num.stems, p=actual.p, simulate.p.value=T, B=9999)
out.sim